R Notebook

Setup

Run AlphaDiversity in scratchnotebooks

#source(here::here("RScripts", "InitialProcessing_3.R"))

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

nonSpikes %>% head

raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

seqDep <- countMat %>% apply(1, sum)
min(seqDep)

[1] 852

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2  3-1-S-180   3-1-S-20    3-1-S-5 
  4.444216   5.154575   4.645763   5.915606   5.132178   3.869332   5.560001   4.583410   4.902727   4.813431   5.336298   4.818610 
 3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20    3-2-B-5  3-2-B-500   3-2-B-53  3-2-S-0-2  3-2-S-1-2  3-2-S-180   3-2-S-20    3-2-S-5 
  4.344849   4.712513   4.675627   5.098950   5.442812   5.030542   4.890340   3.794814   4.924798   4.675124   4.995207   5.202696 
 3-2-S-500   3-2-S-53  3-3-B-0-2  3-3-B-1-2  3-3-B-180   3-3-B-20    3-3-B-5  3-3-B-500   3-3-B-53  3-3-S-180   3-3-S-20  3-3-S-500 
  4.883313   4.402086   4.437311   4.927824   3.459281   5.690676   5.258778   5.316690   5.477329   5.142843   4.917158   4.903164 
  3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180   4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53  4-3-O-1-2  4-3-O-180    4-3-O-5  4-3-O-500 
  4.361907   4.366410   4.965965   4.564515   4.452627   4.567276   4.231160   4.336562   5.019498   4.691336   5.147143   3.990599 
  4-3-O-53  4-3-S-0-2  4-3-S-180   4-3-S-20  4-3-S-500   4-3-S-53  5-1-S-1-2  5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53 
  4.550053   2.805419   4.479398   4.760915   4.733578   4.575868   4.459307   4.512310   4.168196   4.078014   4.011818   4.117664 
 5-5-B-0-2  5-5-B-180    5-5-B-5  5-5-B-500   5-5-B-53  5-5-S-180    5-5-S-5  5-5-S-500   5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 
  4.655036   5.151959   5.481992   5.093104   4.901198   4.294930   4.967999   4.952305   4.203941   4.368584   4.829209   4.897388 
 C_5P1B_20 C_5P1B_500  C_5P1B_53 
  5.468133   4.727174   5.189507 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2   3-1-S-180    3-1-S-20 
0.010656840 0.006422546 0.011469251 0.002865951 0.006273652 0.041830621 0.007109976 0.008699531 0.007976329 0.008690465 0.007414011 
    3-1-S-5   3-2-B-0-2   3-2-B-1-2   3-2-B-180    3-2-B-20     3-2-B-5   3-2-B-500    3-2-B-53   3-2-S-0-2   3-2-S-1-2   3-2-S-180 
0.009319248 0.011267761 0.009263164 0.012319590 0.006327012 0.007362969 0.008080497 0.010290423 0.020464414 0.008360626 0.010250272 
   3-2-S-20     3-2-S-5   3-2-S-500    3-2-S-53   3-3-B-0-2   3-3-B-1-2   3-3-B-180    3-3-B-20     3-3-B-5   3-3-B-500    3-3-B-53 
0.007953267 0.008159670 0.008997352 0.015550953 0.012446876 0.007485125 0.061772871 0.005399056 0.008242407 0.006453177 0.007146198 
  3-3-S-180    3-3-S-20   3-3-S-500    3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180    4-3-B-20     4-3-B-5   4-3-B-500    4-3-B-53 
0.008102995 0.010257642 0.010067304 0.013101847 0.010688464 0.009291003 0.009253490 0.009804310 0.008549470 0.010463207 0.009012522 
  4-3-O-1-2   4-3-O-180     4-3-O-5   4-3-O-500    4-3-O-53   4-3-S-0-2   4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53   5-1-S-1-2 
0.009845773 0.012908904 0.008116955 0.012923854 0.008953551 0.124685269 0.011208241 0.008635149 0.009226702 0.012415445 0.009630661 
  5-1-S-180    5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53   5-5-B-0-2   5-5-B-180     5-5-B-5   5-5-B-500    5-5-B-53   5-5-S-180 
0.012650427 0.013618132 0.015120557 0.012654152 0.012996712 0.008325702 0.008868699 0.007205562 0.009797355 0.009640319 0.013719828 
    5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2  C_5P1B_180  C_5P1B_1P2   C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.009890166 0.010941140 0.016917268 0.007398312 0.008438195 0.005447460 0.003736643 0.006518782 0.005323390 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-226.81  -37.14    4.33   44.82  206.47 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          283088.268 134955.453   2.098 0.039602 *  
log(Size_Class)          30.953      8.177   3.786 0.000324 ***
I(log(Size_Class)^2)     -6.124      1.521  -4.026 0.000144 ***
lat                  -14744.824   7032.354  -2.097 0.039687 *  
I(lat^2)                192.061     91.548   2.098 0.039576 *  
depth                     4.929      3.238   1.522 0.132610    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 76.87 on 69 degrees of freedom
Multiple R-squared:  0.2527,    Adjusted R-squared:  0.1986 
F-statistic: 4.667 on 5 and 69 DF,  p-value: 0.0009955

Richness

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-546.23 -122.72    1.00   96.67 1310.85 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)           1.031e+06  4.685e+05   2.202  0.03105 * 
log(Size_Class)       7.652e+01  2.838e+01   2.696  0.00881 **
I(log(Size_Class)^2) -1.637e+01  5.281e+00  -3.101  0.00279 **
lat                  -5.378e+04  2.441e+04  -2.203  0.03095 * 
I(lat^2)              7.008e+02  3.178e+02   2.205  0.03076 * 
depth                 2.267e+01  1.124e+01   2.016  0.04765 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 266.8 on 69 degrees of freedom
Multiple R-squared:  0.1903,    Adjusted R-squared:  0.1316 
F-statistic: 3.244 on 5 and 69 DF,  p-value: 0.01091

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-464.71 -138.96  -20.87  111.98 1409.67 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           574.310     48.667  11.801  < 2e-16 ***
log(Size_Class)        77.285     28.902   2.674  0.00927 ** 
I(log(Size_Class)^2)  -16.870      5.376  -3.138  0.00247 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 272 on 72 degrees of freedom
Multiple R-squared:  0.1218,    Adjusted R-squared:  0.09738 
F-statistic: 4.992 on 2 and 72 DF,  p-value: 0.009327

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3356 -0.1482  0.0374  0.3250  0.7495 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.623e+03  7.861e+02   2.065   0.0427 *  
log(Size_Class)       1.998e-01  4.763e-02   4.195 7.96e-05 ***
I(log(Size_Class)^2) -3.705e-02  8.861e-03  -4.182 8.36e-05 ***
lat                  -8.432e+01  4.096e+01  -2.059   0.0433 *  
I(lat^2)              1.098e+00  5.332e-01   2.059   0.0433 *  
depth                 1.897e-02  1.886e-02   1.006   0.3180    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4477 on 69 degrees of freedom
Multiple R-squared:  0.2971,    Adjusted R-squared:  0.2461 
F-statistic: 5.832 on 5 and 69 DF,  p-value: 0.0001503

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.014133 -0.004666 -0.002454  0.000611  0.100308 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.507e+01  2.635e+01  -0.952   0.3447  
log(Size_Class)      -3.775e-03  1.596e-03  -2.365   0.0209 *
I(log(Size_Class)^2)  6.540e-04  2.970e-04   2.202   0.0310 *
lat                   1.306e+00  1.373e+00   0.951   0.3449  
I(lat^2)             -1.699e-02  1.787e-02  -0.950   0.3452  
depth                -3.932e-04  6.322e-04  -0.622   0.5361  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01501 on 69 degrees of freedom
Multiple R-squared:  0.09128,   Adjusted R-squared:  0.02544 
F-statistic: 1.386 on 5 and 69 DF,  p-value: 0.2402

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14 
235.0858 235.0858 202.6115 202.6115 216.5492 165.3419 165.3419 196.1950 206.0323 306.2069 306.2069 273.7326 273.7326 287.6703 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28 
236.4630 236.4630 267.3161 267.3161 334.7205 334.7205 302.2462 302.2462 316.1839 264.9766 264.9766 295.8297 305.6670 305.6670 
      29       30       31       32       33       34       35       36       37       38       39       40       41       42 
338.5321 338.5321 306.0579 306.0579 319.9956 319.9956 268.7882 268.7882 299.6413 299.6413 327.1209 294.6466 294.6466 308.5843 
      43       44       45       46       47       48       49       50       51       52       53       54       55       56 
308.5843 257.3769 257.3769 257.3769 288.2300 288.2300 298.0673 298.0673 296.3515 296.3515 263.8772 263.8772 277.8150 277.8150 
      57       58       59       60       61       62       63       64       65       66       67       68       69       70 
226.6076 226.6076 226.6076 257.4607 257.4607 267.2980 267.2980 256.5983 224.1240 224.1240 238.0617 238.0617 186.8544 186.8544 
      71       72       73       74       75 
186.8544 217.7075 217.7075 227.5448 227.5448 

$se.fit
 [1] 26.89171 26.89171 26.15853 26.15853 26.01893 27.68111 27.68111 29.23541 33.91547 18.86417 18.86417 18.25002 18.25002 17.13889
[15] 19.96813 19.96813 21.42573 21.42573 19.18476 19.18476 18.72459 18.72459 17.17411 20.01466 20.01466 21.12886 28.07566 28.07566
[29] 19.18259 19.18259 18.69162 18.69162 16.90082 16.90082 19.54040 19.54040 20.50140 20.50140 18.49526 17.85851 17.85851 15.94069
[43] 15.94069 18.37681 18.37681 18.37681 19.36426 19.36426 26.45911 26.45911 19.24276 19.24276 18.35897 18.35897 16.62051 16.62051
[57] 18.36048 18.36048 18.36048 19.42883 19.42883 26.02108 26.02108 24.37188 23.42219 23.42219 22.26479 22.26479 23.05780 23.05780
[71] 23.05780 24.05365 24.05365 29.12182 29.12182

$df
[1] 69

$residual.scale
[1] 76.86895

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14 
505.4108 505.4108 379.2309 379.2309 480.8810 295.0839 295.0839 434.4257 374.1099 684.3823 684.3823 558.2024 558.2024 659.8525 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28 
474.0554 474.0554 613.3972 613.3972 751.7117 751.7117 625.5318 625.5318 727.1819 541.3848 541.3848 680.7266 620.4108 620.4108 
      29       30       31       32       33       34       35       36       37       38       39       40       41       42 
753.2514 753.2514 627.0715 627.0715 728.7216 728.7216 542.9245 542.9245 682.2663 682.2663 716.6596 590.4797 590.4797 692.1298 
      43       44       45       46       47       48       49       50       51       52       53       54       55       56 
692.1298 506.3326 506.3326 506.3326 645.6744 645.6744 585.3586 585.3586 626.7627 626.7627 500.5828 500.5828 602.2329 602.2329 
      57       58       59       60       61       62       63       64       65       66       67       68       69       70 
416.4358 416.4358 416.4358 555.7776 555.7776 495.4618 495.4618 514.1003 387.9204 387.9204 489.5705 489.5705 303.7734 303.7734 
      71       72       73       74       75 
303.7734 443.1152 443.1152 382.7994 382.7994 

$se.fit
 [1]  93.34957  93.34957  90.80446  90.80446  90.31987  96.08980  96.08980 101.48527 117.73123  65.48346  65.48346  63.35153
[13]  63.35153  59.49448  69.31563  69.31563  74.37542  74.37542  66.59631  66.59631  64.99893  64.99893  59.61672  69.47716
[25]  69.47716  73.34491  97.45941  97.45941  66.58878  66.58878  64.88447  64.88447  58.66804  58.66804  67.83086  67.83086
[37]  71.16680  71.16680  64.20284  61.99249  61.99249  55.33512  55.33512  63.79168  63.79168  63.79168  67.21942  67.21942
[49]  91.84787  91.84787  66.79767  66.79767  63.72976  63.72976  57.69500  57.69500  63.73497  63.73497  63.73497  67.44355
[61]  67.44355  90.32734  90.32734  84.60242  81.30576  81.30576  77.28806  77.28806  80.04086  80.04086  80.04086  83.49774
[73]  83.49774 101.09096 101.09096

$df
[1] 69

$residual.scale
[1] 266.8362

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14 
4.531951 4.531951 4.360594 4.360594 4.338910 4.039890 4.039890 4.161410 4.400479 4.984722 4.984722 4.813365 4.813365 4.791681 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28 
4.492661 4.492661 4.614181 4.614181 5.175134 5.175134 5.003777 5.003777 4.982092 4.683073 4.683073 4.804592 5.043662 5.043662 
      29       30       31       32       33       34       35       36       37       38       39       40       41       42 
5.215580 5.215580 5.044223 5.044223 5.022539 5.022539 4.723519 4.723519 4.845039 4.845039 5.158760 4.987403 4.987403 4.965719 
      43       44       45       46       47       48       49       50       51       52       53       54       55       56 
4.965719 4.666699 4.666699 4.666699 4.788219 4.788219 5.027288 5.027288 4.987930 4.987930 4.816573 4.816573 4.794889 4.794889 
      57       58       59       60       61       62       63       64       65       66       67       68       69       70 
4.495869 4.495869 4.495869 4.617389 4.617389 4.856458 4.856458 4.760224 4.588867 4.588867 4.567183 4.567183 4.268163 4.268163 
      71       72       73       74       75 
4.268163 4.389683 4.389683 4.628752 4.628752 

$se.fit
 [1] 0.15663934 0.15663934 0.15236870 0.15236870 0.15155556 0.16123742 0.16123742 0.17029094 0.19755145 0.10988038 0.10988038
[12] 0.10630304 0.10630304 0.09983095 0.11631071 0.11631071 0.12480098 0.12480098 0.11174773 0.11174773 0.10906736 0.10906736
[23] 0.10003607 0.11658176 0.11658176 0.12307179 0.16353561 0.16353561 0.11173510 0.11173510 0.10887529 0.10887529 0.09844420
[34] 0.09844420 0.11381929 0.11381929 0.11941696 0.11941696 0.10773152 0.10402259 0.10402259 0.09285161 0.09285161 0.10704160
[45] 0.10704160 0.10704160 0.11279330 0.11279330 0.15411952 0.15411952 0.11208560 0.11208560 0.10693769 0.10693769 0.09681145
[56] 0.09681145 0.10694644 0.10694644 0.10694644 0.11316940 0.11316940 0.15156809 0.15156809 0.14196175 0.13643000 0.13643000
[67] 0.12968835 0.12968835 0.13430751 0.13430751 0.13430751 0.14010812 0.14010812 0.16962930 0.16962930

$df
[1] 69

$residual.scale
[1] 0.4477477

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9          10          11 
0.018866864 0.018866864 0.021537876 0.021537876 0.020512033 0.024377201 0.024377201 0.021531485 0.018827124 0.010430513 0.010430513 
         12          13          14          15          16          17          18          19          20          21          22 
0.013101525 0.013101525 0.012075681 0.015940850 0.015940850 0.013095134 0.013095134 0.006715502 0.006715502 0.009386514 0.009386514 
         23          24          25          26          27          28          29          30          31          32          33 
0.008360671 0.012225839 0.012225839 0.009380123 0.006675762 0.006675762 0.005657637 0.005657637 0.008328649 0.008328649 0.007302806 
         34          35          36          37          38          39          40          41          42          43          44 
0.007302806 0.011167974 0.011167974 0.008322258 0.008322258 0.006418762 0.009089774 0.009089774 0.008063930 0.008063930 0.011929098 
         45          46          47          48          49          50          51          52          53          54          55 
0.011929098 0.011929098 0.009083382 0.009083382 0.006379021 0.006379021 0.009130701 0.009130701 0.011801713 0.011801713 0.010775870 
         56          57          58          59          60          61          62          63          64          65          66 
0.010775870 0.014641038 0.014641038 0.014641038 0.011795322 0.011795322 0.009090960 0.009090960 0.012896424 0.015567436 0.015567436 
         67          68          69          70          71          72          73          74          75 
0.014541593 0.014541593 0.018406761 0.018406761 0.018406761 0.015561045 0.015561045 0.012856683 0.012856683 

$se.fit
 [1] 0.005250060 0.005250060 0.005106922 0.005106922 0.005079668 0.005404174 0.005404174 0.005707619 0.006621306 0.003682846
[11] 0.003682846 0.003562945 0.003562945 0.003346021 0.003898371 0.003898371 0.004182938 0.004182938 0.003745434 0.003745434
[21] 0.003655596 0.003655596 0.003352896 0.003907456 0.003907456 0.004124981 0.005481202 0.005481202 0.003745011 0.003745011
[31] 0.003649159 0.003649159 0.003299541 0.003299541 0.003814866 0.003814866 0.004002483 0.004002483 0.003610823 0.003486512
[41] 0.003486512 0.003112095 0.003112095 0.003587699 0.003587699 0.003587699 0.003780478 0.003780478 0.005165604 0.005165604
[51] 0.003756759 0.003756759 0.003584217 0.003584217 0.003244817 0.003244817 0.003584510 0.003584510 0.003584510 0.003793084
[61] 0.003793084 0.005080088 0.005080088 0.004758113 0.004572707 0.004572707 0.004346748 0.004346748 0.004501567 0.004501567
[71] 0.004501567 0.004695985 0.004695985 0.005685443 0.005685443

$df
[1] 69

$residual.scale
[1] 0.0150071

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013938 -0.005151 -0.002547  0.000861  0.105939 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.962e+01  2.755e+01  -0.712   0.4788  
log(Size_Class)      -3.288e-03  1.669e-03  -1.970   0.0529 .
I(log(Size_Class)^2)  5.742e-04  3.106e-04   1.849   0.0688 .
lat                   1.020e+00  1.436e+00   0.711   0.4797  
I(lat^2)             -1.325e-02  1.869e-02  -0.709   0.4806  
depth                -2.376e-04  6.612e-04  -0.359   0.7204  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01569 on 69 degrees of freedom
Multiple R-squared:  0.06714,   Adjusted R-squared:  -0.0004537 
F-statistic: 0.9933 on 5 and 69 DF,  p-value: 0.4284

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7             8             9 
 0.0117020065  0.0117020065  0.0135963432  0.0135963432  0.0133971262  0.0160361854  0.0160361854  0.0140005291  0.0099098120 
           10            11            12            13            14            15            16            17            18 
 0.0043418030  0.0043418030  0.0062361398  0.0062361398  0.0060369228  0.0086759820  0.0086759820  0.0066403257  0.0066403257 
           19            20            21            22            23            24            25            26            27 
 0.0011173144  0.0011173144  0.0030116511  0.0030116511  0.0028124341  0.0054514933  0.0054514933  0.0034158370 -0.0006748801 
           28            29            30            31            32            33            34            35            36 
-0.0006748801  0.0002246683  0.0002246683  0.0021190051  0.0021190051  0.0019197880  0.0019197880  0.0045588472  0.0045588472 
           37            38            39            40            41            42            43            44            45 
 0.0025231909  0.0025231909  0.0009183182  0.0028126550  0.0028126550  0.0026134379  0.0026134379  0.0052524972  0.0052524972 
           46            47            48            49            50            51            52            53            54 
 0.0052524972  0.0032168409  0.0032168409 -0.0008738762 -0.0008738762  0.0033312013  0.0033312013  0.0052255380  0.0052255380 
           55            56            57            58            59            60            61            62            63 
 0.0050263210  0.0050263210  0.0076653802  0.0076653802  0.0076653802  0.0056297239  0.0056297239  0.0015390069  0.0015390069 
           64            65            66            67            68            69            70            71            72 
 0.0066640405  0.0085583773  0.0085583773  0.0083591603  0.0083591603  0.0109982195  0.0109982195  0.0109982195  0.0089625632 
           73            74            75 
 0.0089625632  0.0048718461  0.0048718461 

$se.fit
 [1] 0.005490363 0.005490363 0.005340673 0.005340673 0.005312171 0.005651530 0.005651530 0.005968865 0.006924373 0.003851415
[11] 0.003851415 0.003726026 0.003726026 0.003499173 0.004076805 0.004076805 0.004374397 0.004374397 0.003916868 0.003916868
[21] 0.003822918 0.003822918 0.003506363 0.004086305 0.004086305 0.004313787 0.005732084 0.005732084 0.003916425 0.003916425
[31] 0.003816186 0.003816186 0.003450566 0.003450566 0.003989478 0.003989478 0.004185682 0.004185682 0.003776096 0.003646094
[41] 0.003646094 0.003254540 0.003254540 0.003751913 0.003751913 0.003751913 0.003953516 0.003953516 0.005402041 0.005402041
[51] 0.003928711 0.003928711 0.003748271 0.003748271 0.003393337 0.003393337 0.003748578 0.003748578 0.003748578 0.003966699
[61] 0.003966699 0.005312611 0.005312611 0.004975899 0.004782006 0.004782006 0.004545704 0.004545704 0.004707610 0.004707610
[71] 0.004707610 0.004910927 0.004910927 0.005945674 0.005945674

$df
[1] 69

$residual.scale
[1] 0.015694

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)